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I . INTRODUCTION 

Multicommodity network flow problems emerge when several 
distinct commodities flow through a common capacitated 
network and share one or more arcs that are subject to joint 
capacity constraints. The objective is usually to find the 
minimum cost flow given demands and supplies of the com- 
modities . Problems of this type are also referred to as 
"multicommodity capacitated transshipment problems" and 
"multicommodity flow problems" (MCFPs) . 

The problem of optimizing MCFPs arises frequently in 
logistic systems. As long as only a single commodity is 
involved, even large-scale problems can now be solved 
routinely by specialized network codes that exploit the pure 
network structure of the problem (e.g., GNET by Bradley, 
Brown, and Graves [Ref. 1]) . Such solvers are not directly 
applicable to the general MCFP, and general linear program 
solvers are usually inapplicable as a result of the large 
constraint matrix encountered with typical MCFPs. 

Because of the importance of the MCFP, much effort has 
been devoted to finding efficient, specialized solution 
techniques. Earlier surveys are given by Kennington [Ref. 2] 
and Assad [Ref. 3] . New and effective decomposition methods 
were recently developed by Staniec [ Ref. 4] and show 
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encouraging results. This paper extends that research by 
deriving and implementing a decomposition for the MCFP based 
on a logarithmic barrier function. 

A. STATEMENT OF THE PROBLEM 

In order to formulate the MCFP as a mathematical program 
the following notation is used : 

G = {I,J} is a network with set of nodes I and set 
of arcs J. 

P is the set of commodities (products) flowing on G. 
i e I is a node of G. 

j e J is an arc of G. 

p € P is a commodity flowing through G. 

Np is an |I| X |J| node-arc incidence matrix for each 
product (Nj^ "^2 = ... ~N|p|). 

N is an |I|*|P| x |J|'|PI matrix with matrices Np 
along the diagonal. Os elsewhere. 

A is a |J| X |J|'|P| matrix (I,I,...,I). 
c = (Cj^, . . . , c I p I ) is a vector of arc costs, length 

Ur |Pl . 

X = (x^ , . . . , X I p I ) is a vector of arc flows, length 

|J| • 1 P| . 

bj^ is a vector of joint capacities with length |J|. 
bj^ is the vector (b^,...,bj^) with length |J|'|P|. 
b2 is the vector of supplies and demands for each 
commodity with length | J| ' |P| • 
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Then the MCFP may be stated as follows : 



(P) 


min cx 


(duals) 


(1) 




s . t . Ax ^ b^ 


(Uj^) 


(2) 




Nx = b 2 


(U 2 ) 


(3) 




0 < X < b^ 


(U 3 ) 


(4) 



For the constraint sets (3) and (4) the abbreviated notation 
X € F will be used. 

The stipulation that binds the flow of several com- 
modities to joint capacities is given by (2) and constitutes 
the set of "complicating constraints". Without those con- 
straints, the problem would reduce to independent, bounded, 
single-commodity flow problems. The duals u^^ corresponding 
to (2) are nonpositive. The constraints (4) are redundant, 
but they ensure that x is bounded when the constraints in (2) 
are relaxed. 

B. SOLUTION METHODOLOGY 

The two basic approaches to solving the MCFP (other than 
a standard primal linear programming method) can be charac- 
terized as either decomposition or partitioning techniques. 
The latter employ a special basis factorization within a 
simplex algorithm in such a way that portions of each 
generated basis maintain characteristics of the pure network 
flow problem. Those method are not investigated here. For 
further detail see Kennington and Helgason [Ref. 5] . 
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Decomposition methods solve MCFP by using a master 
problem that coordinates the solution of single network flow 
subproblems. These methods are attractive, since they may 
require the internal storage of only one commodity at a time. 
This approach can further be divided into resource-directive 
and price-directive algorithms. Both will be stated here for 
later reference. 

The resource-directive method uses a master problem that 
distributes improving capacity allocations to the individual 
commodity subproblems. For this purpose MCFP may be written 
as 



(P') 


min 

y 


min 

X 


cx 


(5) 




s . t . 




Ay = 


^1 


(6) 






X 


- y ^ 


0 


(7) 








Nx = 


b2 


(8) 






0 


VI 

>1 

VI 


^1 


(9) 






0 


VI 

VI 


^1- 


(10) 




Any vector y 


satisfying 


constraints 


(6) and (9) may 


be 


interpreted as 


a "capacity 


allocation" 


which apportions 


the 


capacity of an 


arc across 


the individual commodities . 


The 



inner minimization for a fixed y amounts to a restriction of 
(P) . If its solution is feasible in (P) , it yields an upper 
bound on the optimal solution. 
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A 

The subproblems for a particular allocation vector y are 



of the 


form 








(SPl (J) ) 


V(J) 


= min cx 


(11) 






s . t . 


Nx = b 2 


(12) 








0 < X < y 


(13) 


which 


is a 


set of 


single- commodity 


minimum cost 



problems, solved independently for each commodity. 

The master problem than becomes 

(MPl) min V(y) (14) 

s . t . Ay ^ b^ (15) 

0 ^ y < bj^ . (16) 

The objective function in (14) is piecewise linear and 
convex. In theory, the master problem can be solved by 
subgradient optimization or by a cutting plane algorithm. 

Penalty and barrier function decomposition are examples 
of price-directive decomposition. Before describing them, it 
is worthwhile looking at the Lagrangian relaxation of (P) . 
If the joint capacity constraints are placed into the 
objective function with multipliers u^ ^ 0, the Lagrangian 
dual of (P) is 

(LR) max min L(u.,,x) = cx + (b, - Ax) (17) 

Uj^^O x>0 

St . X e F . 
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This problem corresponds to solving max LR(uj^) where 
LR(u-|^) is defined by 

(LR(Uj^) ) 



min cx - (Ax - b^^) 
xe F 

= min (c - u^A)x + 
xe F 



(18) 



Note that the evaluation of LR(Uj^) requires the solution of 
|P| independent single commodity problems. Furthermore, for 
any fixed u^^ < 0, LR(u^) yields a lower bound on (P) and this 

'k 

bound will be tight if u^ is optimal in (P) : 



LR(u2 ) = cx . 

LR(uj^) is a piecewise linear and concave function that 
can be optimized, in theory, by subgradient optimization (for 
example, see Fisher [Ref. 6], Coffin [Ref. 7] or Sandi [Ref. 
8]), or by a cutting plane method ( see Kelly [Ref. 9]). 
Optimization of LR(Uj^) by a cutting plane algorithm is 
essentially equivalent to solving the MCFP by Dantzig-Wolfe 
decomposition (see Staniec, [Ref. 4]). 

However, even if (18) is solved optimally, it is possible 

* 

that LR(uj^) = cx for u^^ . Furthermore, we may not be 

able to obtain a corresponding primal solution that is 
feasible to problem (P) . 

Solving LR(u^) for any Uj^ generally allows some 
constraints to be violated with penalty -Uj^ (Ax - bj^) . Other 
penalty functions are possible which will, at least in a 
limiting sense, yield optimal primal solutions. The 
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objective function in (P) can be transformed into a nonlinear 
auxiliary function Q(h,x) that typically includes a polynom- 
ial penalty term of the form 



(PP(h)) min Q(h,x) 
x€ F 



cx + I II (Ax - b;L) + l I t 



= cx + q(h,x) 

where | | ' | |^ indicates the t^^ norm, and (Ax 



b]^)^ is the 



( 20 ) 



vector max(0,Ax - b^), i.e. the vector of violations of the 
constraint set (2) . The value of h constitutes a positive, 
increasing sequence of penalty parameters. Furthermore, 
t > 1 is required for this method to ensure convergence. 

The penalty function has some attractive properties, as 
given in Luenberger [ Ref. 10] : Let { h } and { x } be 

sequences of penalty parameters and optimal solutions to 
(PP (h^) ) , respectively, with > h'^ and h^ > 0 . Then 



The algorithm converges to the optimal solution. Thus, 

lim Q(h^,x^) = cx , and x^ — > x 
h^^oo 

Since a feasible solution is usually not obtained until final 
convergence occurs, the penalty approach is classified as an 
exterior point algorithm. 




( 21 ) 




( 22 ) 




(23) 
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In order to obtain intermediate feasible solutions a 
suitable scaling or projection procedure may be used. If we 
model the MCFP to include additional bypass arcs which 
satisfy any undeliverable demand at a high cost from a 
supersource, we can assume that any allocation of capacity 
satisfying (6) and (9) is also feasible in F. One pos- 
sibility to obtain such capacity allocations y from a given 

/\ A A 

vector X € F is given by defining y = CA(x) as 



A 



y 



pj 



>=pjbij 



(Ax) j 



if (Ax-b2)j ^ 0 



(24) 



I^Xpj + (bj^ - Ax) j /|P| if (Ax-bi)j < 0. 

A 

Then, a solution of (SPl (y) ) as defined in (24) yields a 

A 

valid upper bound at the current iterate x. This allocation 

procedure was successfully applied by Staniec [Ref. 4]. 

The idea behind the general barrier function approach is 

just the opposite of the exterior point penalty methods. 

Starting with a feasible solution which lies within the 

interior of the constraint region, a modified objective 

function establishes a barrier against leaving the feasible 

region. For the MCFP, the auxiliary function with respect to 

the joint capacity constraints then becomes 

(BP (h) ) min B(h,x) = cx + b(h,x). (25) 

x€F 
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An ideal barrier term b(x) would take the value zero for 
all interior points and infinity at the boundary. A sequen- 
tial barrier function b(h,x) omits this discontinuity and may 
take the following forms ; 

-1 (26) 

b (h, X) = h Z • (b]^ - Ax) . 

or 

b(h,x) = - h Zj ln{h^ - Ax)^, (27) 

where (bj^ - Ax) j denotes the component of the vector 

(bj^ - Ax) . 

The first expression (26) is called the "inverse barrier 
function" and the second term (27) the "logarithmic barrier 
function" (see, for example, Fiacco and McCormick [Ref. 11] 
or Bazaraa and Shetty [Ref. 12]) . Both functions approach 
the ideal barrier function b(x) as h — » 0 . Any barrier 
algorithm requires the existence of an interior region, i.e. 
it does not work for equality constraints. 

The logarithmic barrier function was first proposed by 
Frisch 1955 [Ref. 13]. It was then derived together with the 
inverse barrier function from the Kuhn-Tucker conditions for 
optimality by Fiacco and McCormick [Ref. 11]. The logarith- 
mic barrier function has obtained recent attention in linear 
programming due to its fundamental properties. Megiddo [Ref. 
14] investigates the properties of a weighted logarithmic 
barrier function for general linear programs that places the 
nonnegativity constraints on x into the barrier term while 
requiring strictly positive values for x. This approach 
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leads to a smooth path through the interior of the constraint 
region towards the optimal solution and theoretically 
validates its use in linear programming. Gill, Murray, 
Saunders, Tomlin and Wright [Ref. 15] established a close 
relationship to the projective linear programming algorithms 
initiated by Karmarkar [Ref. 16] . 

The transformation of the objective function into a 
penalty or barrier function creates a nonlinear programming 
problem which requires an efficient solution method to make 
such transformation profitable. Staniec [Ref. 4] successful- 
ly applied the method of restricted simplicial decomposition 
(RSD) in order to solve the penalty function decomposition. 
The idea was developed by Hearn, Lawphongpanich and Ventura 
[Ref. 17] . RSD solves any pseudoconvex optimization problem 
with linear constraints by generating extreme points in 
linear subproblems while a master problem optimizes the 
original objective function over a simplex derived from a 
fixed number of retained extreme points plus the last 
iterate, i.e. the last solution to the master problem. The 
implementation of this method for the MCFP will be described 
later in more detail since it proves useful for the barrier 
decomposition as well . 

C . TEST PROBLEMS 

A real-world large-scale MCFP should be used in order to 
examine the efficiency of the proposed algorithms. Staniec 
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developed an appropriate problem that describes the 
transshipment of conventional ammunition from production and 
storage locations to overseas debarkation points and theatre- 
of-war locations via capacitated road, rail, sea and air 
transportation links. The product demands are time-phased 
and the objective is to minimize the weighted deviation from 
on-time deliveries. The network contains backlogging arcs 
with graduated penalties and bypass arcs that satisfy 
undeliverable demand at high cost. For more details see 
Staniec [Ref. 4 and 18] . The same network is used to test 
and compare the algorithms for four, ten, and 100 commodity 
problems. The underlying network contains 3,300 nodes and 
10,400 arcs of which about 10% are subject to non-redundant 
capacity constraints. 
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II. THE CONCEPT OF BARRIER FUNCTION DECOMPOSITION 



The original MCFP (P) constitutes a linear program. It 
is only the size of the problem that makes a primal simplex 
algorithm unattractive and computationally expensive. 
Penalty and barrier functions were initially designed for 
nonlinear programming problems where they proved useful in 
converting a constrained nonlinear optimization problem into 
an unconstrained nonlinear problem. For the special case of 
the MCFP, penalty and barrier functions provide good 
decomposition tools. 

The barrier function decomposition for the MCFP retains 
the basic decomposition idea of the penalty function 
approach. The overlapping joint capacity constraints in (2) 
are placed into the barrier term of the objective function 
in order to enable the successive solutions of independent 
single commodity problems in a sequence of subproblems. The 
selection of the logarithmic or inverse barrier function is 
not arbitrary but has an interesting theoretical derivation. 

A. THE DERIVATION OF THE BARRIER FUNCTION 

Fiacco and McCormick derive the barrier function from 
the Kuhn-Tucker sufficiency conditions for constrained 
minima [Ref. 11]. A good and detailed analysis, including 
implementations and numerical results, is further given by 
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Wright [Ref. 19] . The derivation shows that the use of a 
logarithmic barrier function is not arbitrary since it has a 
very natural origin. 

For the purpose of this analysis it is convenient to 

restate the problem (P) in the following form : 

(P'') Min f (x) = cx (28) 

xeF 

s.t. g(x)=b^-Ax>0 (29) 

where the objective function (28) has gradient Vf(x) = c and 
each constraint in (29) has gradient vg j (x) = - a^, the 

negative of the row in A. We associate the dual vari- 
ables with the constraints in (29) and note that in 

comparison with the duals u^^ of problem (P), now takes 

the opposite sign : Uj^ = - . 

Following the derivation of Fiacco and McCormick, we 
assume that there exist points in the neighborhood of the 
optimal solution to (P'') such that strict inequality holds 
for the constraints in (29) i.e., g(x) > 0. Furthermore, we 
allow a perturbation of magnitude h in the Kuhn-Tucker 

sufficiency conditions for optimality. At some point 

★ ★ 

[x(h),u^(h)l near the optimum (x , u^^ ) the following 
conditions have to hold for h small and for all j € J: 

(bj^ - Ax) j > 0 (primal feasibility) (30) 

^Ij (b^ - Ax)j = h > 0 (31) 

(perturbed complementary slackness) 
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(32) 



^Ij ^ 0 

(dual feasiblity) 
c - Zj (-a^) = 0. (33) 

Rewriting (31) as = h / (bj^ - Ax) j and substituting in 

(33) yields : 

c - h Zj [(b;^ - Ax)j]"^ (-a^) = 0. (34) 

Using the notation gj = (b^^ - Ax) j ^0 for all j € J , 
equation (34) is of the general form 

Vg . [x (h) ] 

^B(h,x) = Vf(x(h)) - h Z . ^ = 0 (35) 

gj[x(h)] 

and simply means that the gradient of the objective function 
B(h,x) = f(x) - h Zj ln(b^ - Ax) j (36) 

vanishes at x(h) . This is the logarithmic barrier function! 
Note that no constraint qualifications are necessary since 
all constraints in the MCFP are linear. 

Fiacco and McCormick show that the second order 
sufficiency conditions are also satisfied for B(h,x) at 
[x (h) , Uq^ (h) ] near the optimum x* and prove the existence of 
x(h) satisfying these conditions. It is also shown that the 
Hessian matrix is positive definite for small h. 

The same authors obtain the inverse barrier function 
from a modification in the derivation above that enforces 
nonnegativity in (32) by introducing a variable t such that 
X = The logarithmic barrier function seems to be the 

more fundamental approach. 
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B. PROPERTIES OF THE LOGARITHMIC BARRIER FUNCTION 

The basic properties of the barrier function are again 
given by Fiacco and McCormick [Ref. 11] as well as by Wright 
[Ref. 19] and can be stated as follows for the MCFP : 

Ic If 

For a decreasing sequence (h } and associated minima [x } 
the following conditions hold : 



B(h^,x^) < B(h^ ^,x^ (37) 

• If 

for sufficiently small h and bounded (b^^ - Ax)j, and 

cx^ ^ cx^"^, (38) 

Zj ln(b 3 ^ - Ax) > < Zj ln(bi - Ax) (39) 

lim h^ Z- ln(b-i - Ax) -4 0, and (40) 

lim B(h^,x^) cx*. (41) 

h^^O 

The following small example will illustrate these 



properties. The problem consists of a nonlinear objective 

2 

function and a single constraint : Min 2x s.t. x > 1 and 

★ 

has the optimal solution x =1. The barrier function 

2 

B(h,x) = 2 X - h In(x-l) has a closed form solution 
x(h) = 0.5 + (0.25 + 0.25 h) . 

The solutions are shown in Figure 1 for linearly decreasing 
values of h from 20 to 0.1. The change in x as a function 
of h and the corresponding logarithmic term are depicted in 
Figure 2. 
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Figure 1 : Trajectory of a Barrier Function Optimization 




Figure 2 : x(h) and b(h,x) for h decreasing 
It is interesting to observe that the path of x(h) 
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towards its boundary x = 1 is smooth and can be well 
approximated by a straight line until B(h,x) reaches its 
maximum. From that point on the rate of change in x (h) 
increases slightly. Similar results are obtained for 
different objective functions. Returning to the multiple 
constraint barrier function for the MCFP, the maximum of 
B(h,x) is obtained for a value h' which satisfies the sta- 
tionarity condition Zjln(bj^ - Ax(h'))j= 0. This can only 
occur if some components of (b^ - Ax(h)) are less than or 
equal to one, i.e., some constraints are almost tight. 
Thus, we will not be able to observe the property in (37) 
until the final stage of the algorithm when relatively small 
values of h are attained. 

The properties (38) and (41) are the most important ones 

for practical purposes. Starting with a feasible (interior) 

point, the algorithm produces a nonincreasing sequence of 

objective values which converges to the optimal value of (P) 

and provides intermediate feasible solutions along its path. 

This path of x as a function of h describes a trajectory 

that has already been studied by Fiacco and McCormick [Ref. 

11] as well as by Wright [Ref. 19] . The existence of this 

trajectory could be utilized in an extrapolation technique 

★ 

predicting the next iterate or even x . Its implementation 
for the large-scale MCFP is not investigated here. However, 
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we will be able to get a good feasible solution after only a 
few iterations without using an extrapolation technique. 

One final property is the robustness of the barrier 
function method: Mifflin [Ref. 20] showed that it is 
sufficient to solve B(h,x) at each iteration only 
approximately, within a predetermined tolerance, while still 
achieving convergence (at a lower rate) . 

C. COMPUTATIONAL CONSIDERATIONS 

Any barrier function technique requires an interior 
starting point. It may not be easy to find a starting point 
and the performance of an algorithm is influenced by the 
quality of this initial solution. We utilize the capacity 

A A 

allocation mechanism y = CA(x) for the generation of an 
initial starting point. 

The transformation of (P) into a nonlinear programming 
problem requires an effective NLP solution methodology. If 
a line search is part of this method, any discrete-step line 
search procedure along an improving direction may cause the 
evaluation of B(h,x) at one or more infeasible points. This 
requires extra precautions during the implementation, 
resulting in additional computation time. 

Ryan [Ref. 21] points out that for small values of h the 
auxiliary function B(h,x) becomes very "steep valleyed" and 
the gradient VB(h,x) can take large values in a small 
neighborhood of x(h). Therefore, a termination criterion 
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based on the magnitude of the gradient alone may be critical 
as h becomes small; termination criteria based on the 
difference in successive objective values may lead to 
premature termination. We will consider both criteria and 
take the risk of not solving the problem optimally at each 
step . 

Another well-known difficulty of barrier function 
approaches is the ill-conditioned Hessian matrix of B(h,x) 
for small h (see, for example Bazaraa and Shetty [Ref. 12], 
Wright [Ref. 19] or Ryan [Ref. 21] ) . This problem emerges 
only in the final stage as h — > 0. For large-scale program- 
ming purposes the achievement of a solutiom which defers 
from the optimum only by some small value € ( € -optimality ) 

is normally sufficient. Such a solution may be obtained 
before the ill-conditioning becomes bothersome. 

Another issue that significantly influences the perfor- 
mance of the algorithm is the choice of the initial barrier 

parameter h and its rate of decrease. Lootsma [Ref. 22] 

★ 

showed that the absolute difference | I x(h) - x | | is of 

the order 0(h) for the logarithmic barrier function. Ryan 

[Ref. 21] uses this linear relationship to propose a 

k 1 — k 0 

generating relation of the form h = 10 h , where 
k = 1,2,3,..., and h® is positive. 

The suitable choice of the initial value h® is a more 
critical issue, since theoretically h® can take any value 
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greater than zero and up to infinity. Intuitively however, 
h^ should depend on the cost and constraint structure of the 
problem. One possible choice would be to interpret the 
parameter h as a scaling factor between the cost vector and 
the set of constraints placed into the objective function. 
This leads to the suggestion that we achieve initial balance 
at the starting point such that cx^ = h^ Zj In (bj^-Ax*^) j . 

Instead of analytically deriving h^, simply taking an 
multiple a (a greater than one) of the maximum cost value in 
connection with a constant rate of decrease has performed 
well in test problems : 

h*^ = a * c__„ , = a * h'^ for k > 1 and 0 < a < 1. 

In general, we recommend choosing h*^ too high rather than 
too low as the algorithm will more quickly adjust to high 
values rather than low values. 



20 



Ill . 



THE LOGARITHMIC BARRIER FUNCTION DECOMPOSITION 



USING RESTRICTED SIMPLICIAL DECOMPOSITION (RSD) 

The basic idea of the barrier function decomposition for 
the MCFP is to place the coupling constraints (2) into the 
logarithmic term of the barrier function as derived in 
Chapter 2 . The resulting formulation (BP (h) ) constitutes a 
nonlinear programming problem with linear network flow 
constraints (3) and (4) . Using the restricted simplicial 
decomposition technique (RSD) , we will decompose the problem 
into a nonlinear master problem and a set of subproblems, 
which require only the solution of |P| independent pure 
network flow problems. The master problem has a reduced 
search space described by a fixed number of retained extreme 
points, which are generated in the subproblems. 

Lower and upper bounds on the optimal solution to (P) 
can be easily established. The analogy with the penalty 
decomposition is interesting and worth pursuing. The 
solution method RSD used in the penalty function decomposi- 
tion also proves to be effective for (BP (h) ) and is descri- 
bed first. 

A. RESTRICTED SIMPLICIAL DECOMPOSITION (RSD) 

The basic difference between a linear programming 
problem and a nonlinear programming problem with linear 
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constraints is the fact that the optimal solution will 
normally not be an extreme point of the constraint region. 

The familiar Frank-Wolfe algorithm [Ref. 23] takes 
advantage of the specialized constraints by generating an 
extreme point solution of the original problem in a linear 
subproblem whose objective function is the linearization of 
the original objective function at the current iterate 
(given an initial solution) . A master problem provides a 
new iterate via a simple line search between the previous 
iterate and the new extreme point. The main disadvantage of 
this decomposition algorithm is its susceptibility to slow 
convergence, especially when the line search direction 
becomes orthogonal to the gradient of the objective function 
(e.g., see Wolfe [Ref. 24]). 

The method of simplicial decomposition is due to 
Geoffrion [Ref. 25], von Hohenbalken [Ref. 26] and Holloway 
[Ref. 27]. A nonlinear master problem replaces the line 
search of the Frank-Wolfe method by extending the optimiza- 
tion to the convex hull of all extreme points generated in 
the linear subproblems. This solution method relies on an 
effective solution of the nonlinear master problem. 
Although the master problems have only simple convexity 
constraints, the increasing number of extreme points makes 
the implementation of this technique unattractive for the 
solution of large-scale programming problems. 
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The disadvantage of simplicial decomposition led to the 
idea of restricted simplicial decomposition (RSD) as develo- 
ped by Hearn, Lawphongpanich and Ventura [Ref. 17] . RSD 
limits the size of the master problem by fixing the maximum 
number r of retained extreme points . The master problem 
optimizes over the simplex of these extreme points and the 
iterate obtained from the previous master solution. Any 
newly generated extreme point replaces the old extreme point 
with minimum weight in the expression of the current iterate 
as a convex combination of the retained extreme points and 
the prior iterate. After solving the new master problem, 
all extreme points with zero weight can be discarded. 

If r is set to its minimum value r = 1, RSD specializes 
to the algorithm of Frank-Wolfe. For the maximum value of r 
(the finite number of extreme points) the method represents 
simplicial decomposition. The solution to the master 
problem becomes harder as r is increased, but is rewarded by 
significant improvements in the convergence rate to the 
optimal solution (see Hearn, et al . , [Ref. 28]). 

The decomposition of (P'') is formed in the following 
k. t h. 

manner: Let X denote a matrix in the k iteration of the 
master problem whose columns are a set of r extreme points 
from F, let x be the solution of the previous master 
problem, and let U {x*^“^}. The master problem in 
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terms of the weights w at iteration k becomes, 
for a fixed h, 

(MP2^) min B(h,X^w) = cX^w - h In (b^j^ - AX*^w)j{42) 



s . t . 



1 w =1 
w > 0, 



(43) 

(44) 



Ai. '^k'^k’^'k 

which has solution w in terms of w and solution x = X w 
in terms of x. 

The subsequent subproblem optimizes the linear ap- 
■ '^k 

proximation of B(h,x) at x over F, which is equivalent to : 

:k. 



(SP2^) 



min vB(h,x )x . 
xeF 



(45) 



It is convenient to introduce the notation (b^ - Ax)^ to 
represent the vector (b^ - Ax) whose components are taken to 
the t power. Then, (SP2 ) can be written as 



(SP2^) 



min (c - h[(bj^ - Ax^) A) x 

xeF 



(46) 



Using the relationship in (31), we observe that we obtain an 

A A _ 1 

estimate H of the optimal dual variables as h(b 2 ~Ax ) 

at each iteration. Substituting into (46) yields a sub- 
problem as 

(47) 



(SP2^) 



min (c + H^A)x. 
xeF 



This subproblem resembles a standard Dantzig-Wolfe decom- 
position subproblem as it would be used in a dynamic column 
generation approach to problem (P) with dual estimates 

A A . 

Ui = - ^^l• This issue is not discussed here; see Staniec 
[Ref. 4] . 
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Due to the block structure of the constraint matrix N, 
k 

(SP2 ) permits independent solutions for each commodity. 
However, each solution to (SP2'^) is not necessarily a 
feasible flow in (P) and of course, is not an interior 
point. As a starting point, the barrier function approach 
requires an initial interior point, which may be obtained as 
follows. First we solve a subproblem with h = 0, and 

yielding solution x® . Then, using the capacity allocation 
mechanism ^ = CA(x^), we further reduce the capacity by 

A A 

setting y' = b * y for 0 < b < 1. Finally, we solve 

A 

subproblem (SPl(y')) to get an initial interior point 
solution to (P) . 

Starting with this solution in the master problem we 
simply limit our search to that part of the convex hull of 
extreme points plus the current iterate which provides an 
interior point solution as next iterate. This will be 
discussed in connection with the solution of the nonlinear 
master problem. 

B. SOLUTION OF THE MASTER PROBLEM 

The reduced gradient method is a well-known approach to 
solving a nonlinear program with linear constraints and is 
used here for solving the RSD master problem. The method 
was first proposed by Wolfe [Ref. 28] and modified by 
McCormick [Ref. 29] . The basic idea is the partition of the 
variables x into m basic variables Xg and n nonbasic 
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variables Xj^ as done in the simplex algorithm of linear 
programming. This induces a partition of the constraint 
matrix A into parts B and C, where B is assumed nonsingular. 
The NLP then takes the form 



min f(x) = f(Xg,Xj^) 




(48) 


s . t . Ax = Bxg + CXj^ = 


b 


(49) 


^B'^N - 


0. 


(50) 



The variables Xj^ are regarded as independent variables 

whereas Xg are dependent variables completely determined by 
Xj^ and equations (49) . Consequently, the objective function 
can be considered to be a function of Xj^ only and the 
constraints reduce to the nonnegativity constraints on the 
independent variables and the limitation in their change 
that provides nonnegative basic variables. This fact allows 
the application of a modified steepest descent method 

accounting only for the nonnegativity constraints. The 

• Ic • • * • 

reduced gradient r at iteration k is of dimension n-m and 

computed as : 

^ ~ ^ ~ ^x ~ ^x ^ ^^B' ^N^ ® 

N B 



To find an improving direction d such that 7f(x^)d < 0, we 
select at each iteration k 



/^ - r. 



dx, • = 

N j ^ 



V_x.r. 



if rj ^ 0 



if rj > 0 



for Xj nonbasic 



(52) 
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and 



dgj^ = - (B~^Cdjj) £ for basic. (53) 

A new direction needs to be computed as soon as a nonbasic 
variable attains its zero level. If a basic variable 

becomes zero, the partition must be modified. Also, this 
method requires a nondegenerate solution at each iteration. 
A more detailed description can be found in Luenberger 
[Ref. 10] or Bazaraa and Shetty [Ref. 12] . 

The use of the reduced gradient method for the solution 

If 

of the master problem (MP2 ) results in some nice simplif- 
ications due to the presence of only a single convexity con- 
straint. There is only one basic variable to be selected 
and the reduced gradient becomes 

rj = cxj - + h[(bj^ - AXw)”^]"^ (AXj - (54) 

for the nonbasic variables w j . Computing the direction 
component for the basic variable reduces to 



di - - B-lcd^ - - Ijd^ 
and for the nonbasic variables Wj 



r - r. 



if rj ^ 0 



(55) 



(56) 



^ “Wjrj if rj > 0. 

The line search in the direction d is limited by the non- 
negativity constraints. Thus a maximum steplength ot' is 



27 



computed as 

a' = min I dj^ < 0 , for all variables Wj^} . (57) 

Any move in the direction d of size a must also satisfy 
- [AX(w + ad)] >0 for all 0 ^ a ^ a' (58) 
or equivalently 

b^ - AXw > (XAXd for all 0 ^ a ^ a' . (59) 

Since b^^ - AXw > 0 , this holds for all arcs 1 with 
(AXd) <0. If (AXd) > 0 for some arcs 1, we perform a 
ratio test to select 

a" < min { (b^^ - AXw) / (AXd) | (AXd) > 0) (60) 

and set 

“max {a' ,a" } . (61) 

The master problem solution is summarized as follows: 
Step 0 : Set w such that Iw =1, Wj^ > 0 . 

Select a basic variable: use the largest Wj^. 

Step 1 : Compute the direction d as determined in equations 
(54) through (56) . If d = 0, stop. 

Step 2 : Solve min B[h,X(w+ad)j s.t. 0 ^ a ^ “max ^ 

line search, yielding Let w <— w + ttj^d and 

go to Step 1 • 

The convergence of this method may not be satisfactory 
since it represents a greedy descent method in the reduced 
space of the nonbasic variables. Reklaitis et al. [Ref. 30] 
suggest a convergence acceleration technique by using a more 
effective unconstrained search method as long as the basic 
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variables do not change. This is especially relevant in our 
case where the prior iterate frequently takes the largest 
weight and remains the basic variable. Therefore we apply a 
conjugate gradient method, that has proven its effectiveness 
in unconstrained optimization and is easily implemented. 
The modified direction for the first iteration and any 
iteration following a basis change becomes : 

- rj if Wj > 0 or rj < 0 



d 



k 

Nj - 



(62) 



0 otherwise, 

and for all other iterations : 



_ _k . 

d N - - r + 



II I 1 2 



k-1 



I l|2 ^ 



(63) 



The direction for the basis variable is the same as in the 
standard reduced gradient method. The use of the conjugate 
gradient provides significant improvements. In a typical 
master problem with only four extreme points, the reduced 
gradient method used 500 iterations to obtain a solution 
that was still 13% worse than a solution obtained with the 
conjugate gradient method in 28 iterations. 

Each line search requires frequent evaluations of the 
objective function and consumes a fair amount of computation 
time. We can improve this by modifying the objective 
function to include only those arcs in the barrier term that 
are potentially able to violate the joint capacity con- 
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straints at the current iteration. These arcs are easily 
identified and recorded in a set containing all arcs that 
ever had a violation in any of the extreme points generated 
in the subproblems. The arcs j i Jv cannot violate the 
joint capacity constraints in a solution to the master 
problem that is a convex combination of the extreme points 
and the feasible prior iterate. Thus we establish the 
barrier only on a reduced subset of the joint capacity 
constraints. This set is updated at each solution of a 

subproblem in case a new arc experiences a capacity viola- 
tion. This procedure amounts to a modified barrier func- 
tion. However, in a finite number of iterations, the number 
of arcs in will take a fixed value IJyl ^ IJI ^nd no more 
arcs will be added to From this point on we are back to 

a pure barrier function as derived in Chapter 2 and conver- 
gence of the algorithm is preserved. 

C. LOWER AND UPPER BOUNDS 

Since we will rarely be able to find the optimal 
solution to (P) within a reasonable number of iterations, we 
need to establish bounds on the optimal solution. 

Lower bounds can be derived from the Lagrangian dual 
problem 

(LR(uj^)) min (c - u-j^A) x + u^bj^ , which provides a lower 
X€ F 

bound on (P) for any fixed Uj^ ^ 0. 



30 



, A A 

Using the dual estimates v.i.z. 

^ A A 

= - n (b^ - Ax) for each solution x provided by the 
master problem, we obtain 

(LR(jx^)) min (c + x - . (64) 

xe F 

Recalling the subproblem (SP2^) : min (c + M-^^A) x , we find 

x€ F 

by comparison that both objective functions differ only by 

A 

the constant term and yield the same optimal solution 

k ^ k 

X . Thus, we obtain a valid lower bound V(|i^^) by subtract- 

A 

ing the constant term M-^b^ : 

V(G^^) = (c + Uj^’^A)x^ - (65) 

k * k * 

Furthermore, since in the limit x — > x and — » |i^ , it 

must be that V(ji^^) — > (c + A)x - |i]^*b||^ = cx* . 

Upper bounds on the optimal solution are generated in 
each master problem, since we restrict the solution to an 
interior, feasible point. Thus, if x*^ = X^w^ solves (MP2^) , 
the upper bound 

V(x^) = cx’^ (66) 

is readily available at each master problem solution. Due 
to the convergence property (41) of the barrier function, 

A T_ ★ 

it must be that V(x ) — > cx as h — » 0. 

However, we are usually able to obtain a better upper 
bound by utilizing the capacity allocation mechanism 

A A k at- 

again. Let y = CA(x ) at some iteration k and let x be an 
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interior point. Then, for all arcs j e Jy, (Ax - b^)j < 0 

A k 

and Ypj > Xpj for all p, j from equations (24). Then 

the following relationship must hold: 



Ay AT- 

V(x^) = CX^ 



min cx 
xeF 
:k 



min cx = V (y ) , 
xeF 



x^x x^y 

A A ^ 

since y ^ x . 

A 

Thus, solving SPl (y) yields a feasible solution with value 
V(y) < V(x*^) . 



D. THE ALGORITHM RSD (B) 

The algorithm RSD (B) using a barrier function decomposi- 
tion can now be presented. An initial lower bound is 
obtained by solving the problem without the joint capacity 
constraints (2). If this solution is feasible in (P) , it is 
optimal. Otherwise we obtain a feasible solution via the 
capacity allocation mechanism, which gives an initial upper 
bound. The algorithm generates a sequence of extreme points 
in the subproblem and an interior point solution in the 
master problem until € -optimality is achieved. As a 
heuristic, we invoke the capacity allocation mechanism at 
every r^^ iteration of the master problem to improve the 
current upper bound and decrement h at this time, even if 
the master problem is not solved optimally for the cur- 
rent h. 



32 



Notation : 

X matrix of retained extreme points at iteration k 

^k 

X current master problem solution 

Xj^ previous master problem solution 

^ Ic • Ic 

X matrix X augmented with x^^ 

)c • V” 

X optimal solution to subproblem (SP2 ) 

H(X) convex hull of X 

CA(x) a capacity allocation based on x 

Jv set of joint capacity arcs which are violated in 
at least one subproblem solution x . 
h^ barrier parameter used in 1^^ parameter update 
r maximal number of retained extreme points 

e stopping criteria for near-optimality. 



Algorithm RSD (B) 



Input : The network G = joint capacity vector b^, 

cost vector c and supply /demand vector \> 2 - 

• • "'k . — 

Output : Best obtained solution x and final bounds V, V. 

Step 0 : (Initialization) 

Select h° > 0, e >0, r > 1, 0 < b < 1. 

Set k = 0, 1 = 0, 

X® = 0 , Jy = 0 , iX^j 0 j € J. 

Solve (LR(0)) : x*^ = argmin {cx | x e F }. 

Set V = cx°. Set Jy = { j I (bj^ -Ax°)j < 0). 
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A Q 

If Jv = 0. stop with X optimal. 

Else, set y = CA(x®), y' = b * y and 

A A 

solve SPl (y' ) yielding x^ . 

Set V = cxy®. If (V - V)/V < e, exit with x^*^, V, V. 
Else, set = 0 , Xj^ = x^ = x^^ , k = 1 . 

Step 1 ; (Solve SP2^) 

Solve x^ = argmin { (c + }i 2 ^^A)x | x e F}, 

X 

where = h^ (b^ - Ax^) for all j 6 Jy and 

^Ij = 0 for all j Jy. 

Set Jy = Jy U {j I (bj^ - Ax^) < 0}. 

Set = (c + - M-i b^. 

If < V, set V = . 

If (V - V) / V < e, exit with V, V. 

If (c + M-i’^A) (x^ - x*^) > 0, x^ solves B(h^,x). 

Set h^"^^ = a * h^ (0 < a ^ 1) . 

Set 1=1+1 and go to Step 2. 

Else 

(i) if |X^| < r , X^'*’^ = x’^ U {x^} . 

(ii) if I X^ I = r , drop the column of X^ which had 

. ^'k . 

the smallest weight w^ in the convex com- 

• • « ^ Ic • • Ic 

bination forming x and replace it with x . 

Go to Step 2. 
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step 2 : (Solve the master problem MP2*^) 

Set U {xj^} . 



A , -l 

Find X = argmin { 


B(h^,x) 


1 X € H{X^'*’^) 


} 


X 










where 


1 ^ i ^ IX^'*’^ 1 


and 


x^ is 


the i^^ 


column in X^'*’^ 


• 



Set Xj^ = . 

If < V , set V = cx^'*’^. 

If (V - V) / V < e, exit with x^, V, V. 

If k+1 is an integer multiple of r, do a capacity 
allocation : Set y = CA(x ) and solve SPl (y) , 
yielding V(y) = cx^^''’^ . 

If V(y) < V, set V = V(y) . 

If (V - V) / V < e, set x^ = , 

exit with x^, V, V. 

Set h^'*'^ = a * h^ (0 < a ^ 1) and 1 = 1 + 1. 

Set k = k + 1 and go to Step 1. 

C. IMPLEMENTATION OF RSD (B) 

The algorithm RSD(B) has been coded in FORTRAN which is 
still an extremely efficient language for mathematical 
programming purposes (see MacLennan [Ref. 31]). A sophis- 
ticated data structure is used for the storage of sparse 
matrices and vectors. Allowing direct communication with 
the X-system solver of Brown and Graves [Ref. 32] . Integer 
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arithmetic is performed to the extent possible, especially 
for the subproblems which use a GNET solver [Ref. 1] to 
produce rapid solutions. 

A design goal was set to provide a decomposition 
procedure which solves only a single commodity problem at a 
time, and thus operates easily within a modest memory 
region (say, two megabyte virtual storage capacity) . Other 
information such as current incumbent, previous iterate, 
prior extreme points, etc., have to be kept on external 
storage devices. This approach leads to considerable 
input/output operations at the expense of computation time, 
but the maximum problem size in terms of number of com- 
modities remains independent of the available virtual 
storage. Also, the resulting algorithm is highly parallel 
by commodities . 

The implementation of the master problem contains 
several parameters such as the number of retained extreme 
points, stopping criteria for optimality at each iteration, 
the final interval of uncertainty in the line search and an 
upper limit on the maximum number of line searches con- 
ducted. Furthermore, the weights of the extreme points and 
the objective function evaluation are subject to roundoff 
errors. For a sensitive objective function, special care 
is necessary to insure convergence. Fortunately, we will 
confirm that the barrier function decomposition is a robust 
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procedure that does not require an optimal solution to the 
master problem at each iteration. Good results are obtained 
over a relative broad range of parameters. 

The comparison of the barrier algorithm RSD(B) with the 
penalty algorithm RSD(P) of Staniec [Ref. 4] reveals that 
they are very similar in their sequential structure. This 
similarity permits embedding both algorithms in the same 
computer program and creates the potential for devising a 
hybrid algorithm which takes advantage of each. The 

relationship between RSD(B) and RSD(P) is easily establish- 
ed. The dual estimates obtained from the penalty approach 
take, for some vector x and joint capacity constraint j, the 
form 

= h [(Ax - b;L)^"^]j , t > 1 (67) 

versus 

= h[(bi - Ax)"^]j , (b^ - Ax)j > 0 ,j e (68) 

for the barrier approach. However, the gradient of the 

penalty function at some vector x takes the same form as for 
the barrier function, namely 

VQ(h,x) = c + jx^A (69) 

versus 

VB(h,x) = c + [i^A. (70) . 

Thus, besides the fact that the barrier function decomposi- 
tion requires an initial interior point, both methods use 
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the same subproblems and master problem 
for the dual estimates and different 
routines in the line search. 



with different input 
function evaluation 
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IV. COMPUTATIONAL RESULTS 



In order to assess the capabilities of the algorithm 
RSD (B) , we solve different versions of the test problem 
described in Chapter I . This problem suite has been exten- 
sively studied by Staniec [Ref. 4], using different al- 
gorithms. The optimal solution for four and ten commodity 
problems is available for comparisons, obtained by solving 
the problem (P) using the X-system [Ref. 32] . The four 
commodity problem (4H) has approximately 13,200 constraints, 
41,600 variables and optimal solution value 130,739,585. The 
ten commodity problem (lOH) has about 33,000 constraints, 
104,000 variables and optimal solution value 169,532,339. 
For purposes of direct comparisons, the penalty algorithm 
RSD(P) of Staniec [Ref. 4] has been converted into our data 
structure and uses our computer program framework. RSD(P) 
has been improved by taking advantage of the conjugate 
gradient modifications in the master problem. 

First, we will evaluate the performance of RSD (B) with 
different initial parameters in solving problem 4H. Subse- 
quently, the comparisons between RSD(P) and RSD (B) are 
presented for problems 4H and lOH. Possible modifications of 
algorithm RSD (B) are then discussed. Finally, we test with 
a 100 commodity problem having more than one million vari- 
ables and 300,000 constraints. 



39 



A. PERFORMANCE OF THE ALGORITHM RSD (B) 

The following results were obtained on an IBM 3033 AP 
under VM/CMS. The Central Processor Unit (CPU) utilization 
time is used as a performance measurement. An "e -optimality 
gap" is computed from the current upper and lower bounds as 
( V - V ) / cx or estimated as ( V - V ) / V . 

Since the barrier function decomposition requires an 
interior starting point, we will first investigate the 
response of RSD(B) to different starting points. Let 
denote the maximum cost value in the network. While fixing 
h® = 2 * '^max' ^ “ 0.5, and r = 7 in problem 4H, the choice 

A A 

of the parameter b establishing y' = b * y resulted in the 
optimality gaps listed in Table 1: 



Gap (%) 



initial gap 


66.76 


51 . 93 


39.42 


26.95 


iteration 7 


10.21 


9.87 


9.56 


10.49 


iteration 16 


4 . 60 


3.98 


4.60 


5.06 


iteration 24 


3.35 


2.25 


3.05 


2.73 


iteration 32 


2.03 


1.42 


1.53 


1.86 


iteration 40 


1.19 


0.99 


1.28 


1.38 


iteration 48 


0.75 


0.68 


0.89 


0.89 


reduction b 


0.80 


0.85 


0.90 


0.95 


Table 1 : 


Response 


to Different 


Starting 


Points 



Problem 4H 
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If the parameter b is selected too low, the dual es- 
timates associated with the resulting interior starting 
point may not generate good, initial extreme points. If the 
starting value is chosen too high, the interior point gets 
too close to the boundary of the constraint region. The sub- 
sequent line search in the master problem is confined to a 
smaller search space of the constraint region resulting in 
slower convergence. Both values, b = 0.85 and b = 0.9, work 
well in the test problem. Further analysis will be based on 
b = 0.9. Good results with RSD were obtained for any value 
of r between 6 and 8 . 

The lower bounds obtained with the algorithm RSD(B) are 

sensitive to the initial parameter h*^ . Recall that the dual 

—1 

estimates are computed as = h(Ax-b 2 ) • If the barrier 
parameter h*^ is relatively small, we approach the boundary 
very soon. As the slack on some jointly capacitated arcs 
almost vanishes, its reciprocal may take huge values result- 
ing in extreme dual estimates. This situation leads to large 
oscillations in subsequent solutions of the subproblem. The 
phenomenum is demonstrated in Figure 3 where h*^ = 0.5 * <^max* 
As h*^ is increased, this effect seems to disappear, as shown 
in Figures 4 through 6. It is interesting to observe that 
some good lower bounds are obtained under all conditions. 

The extreme result obtained for the small value of h® 
suggests that h*^ > is the better choice. 
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Figure 3: V{p.^) for h® - Problem 4H 




Figure 4: ViV^^) for h° = l*c^ax' Problem 4H 
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LOWER BOUND 




Figure 5: V(p,j^) for h® = Problem 4H 




Figure 6; V((lj^) for h® = Problem 4H 
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If we want to select an initial value that balances 

the linear part cx and the barrier term in the objective 
function, we have to take into account that IJyl fixed 

but initially increasing. We do not have sufficient inform- 
ation about the size of the barrier term unless at least one 
pilot run has been conducted. Using the information that 
about 6% of the arcs are finally contained in the set Jy for 
problem 4H, we would obtain a value of h^ < 0.5 c__„, which 

is not the best choice, but may serve as a lower bound on h*^ . 

The differences obtained for the upper bounds are less 

A V 

significant. Figure 7 shows the sequence of values cx 
obtained for the same choices of h*^ as before. 



% 




A V 0 

Figure 7: Response of cx to h , Problem 4H. 
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Any decrease in h creates a disturbance in the dual es- 
timates. Since the master problem is not solved optimally 
for each value of h, the algorithm may not always provide 
better lower bounds before the next update of h occurs. We 
are still able to get good lower bounds when a relatively 
moderate rate of decrease a = 0.5 is used. 

A solution trajectory for problem 4H is given in Figure 8 
using h*^ = 1.5 * where a 0.6% optimal solution is 

obtained within 200 seconds . Instead of plotting the 

A 

current upper and lower bounds, the current values of V()i) 

A V . . 

and cx are shown. We observe the almost strictly decreasing 

]c 

sequence of cx although the master problem is not solved 
optimally. The improved upper bounds obtained by the 
capacity allocation mechanism are denoted as cx' and indi- 
cate that significant gains are obtained only at the initial 
stage. The differences between cx and cx' nearly vanish 
after about 100 seconds although slightly tighter upper 
bounds are produced throughout. 
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Figure 8: Trajectory of RSD (B) , Problem 4H 

It seems to be typical for the barrier function that good 
upper bounds are obtained at an early stage, whereas the 
lower bounds trail behind. 

The trajectory of the objective value B(h,x ) together 

^ Ir 

with its linear part cx is given in Figure 9. We find that 

A ^ at 

B(h,x ) still approaches the optimum cx from below, the 
barrier term takes negative values over the whole range. The 
steps in its trajectory are due to the decreases of h by a 
factor of 0.5. 
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Figure 9: Values of the Objective Function, Problem 4H 

A direct comparison between the algorithms RSD(P) and 
RSD(B) shows that the barrier function decomposition is less 
effective in the solution of the smaller problem 4H but is 
competitive in the problem lOH. The solutions to problem 4H 
are compared in Figure 10 were h*^ « *^ ' *^*^^ **^max RSD(P) and 
h° = RSD(B) . 
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Figure 10: Comparison of RSD(P) to RSD (B) , Problem 4H 

If we use the same initial parameters for the solution of 
problem lOH, we obtain an interesting result. The trajec- 
tories are given in Figure 11. Obviously, the initial 
barrier and penalty parameters, respectively, are too low in 
both cases, yielding poor and oscillating lower bounds for 
RSD(B) versus bad upper bounds for RSD(P) due to insufficient 
penalty on the capacity violations. (This situation suggests 
investigation of a hybrid algorithm, incorporating both 
RSD (B) and RSD (P) ) . 
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Figure 11: RSD(P) Versus RSD (B) , Problem lOH, 

Same Initial Parameters 

The initial parameter h® has been increased until good 
results are obtained for both methods. The result is shown 

. . ''If 

in Figure 12. We observe that the final values of V(p. ) 

decay for both algorithms. We do not generate improving 
lower bounds. Apparently, the parameter h is updated to 
rapidly before the master problem generates good dual 
estimates. Staniec [Ref. 4] reported similiar poor lower 
bounding for problem lOH. 
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Figure 12: RSD(P) Versus RSD(B), Problem lOH, 
Best Parameters 



In summary, we find that RSD(B) provides good upper 
bounds at each iteration and does not depend on the capacity 
allocation routine as the penalty algorithm does. On the 
other hand, the penalty algorithm provides better initial 
dual estimates, resulting in better lower bounds. Both 
algorithms converge to a good solution within a reasonable 
amount of time. 
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B. POSSIBLE MODIFICATIONS OF THE ALGORITHM 

A first modification could be to postpone a further 
decrease of the barrier function parameter h until at least 
one better lower bound has been obtained compared to the 
bound generated for the previous value of h. Note, however, 
that the sequence of upper bounds is almost strictly decreas- 
ing from iteration to iteration and controlled by the 
magnitude of h. A higher barrier parameter results in higher 
upper bounds. Therefore, this choice should depend on the 
problem at hand. If a fast approximation of the solution is 
desired, h should be decreased more rapidly. If a higher 
resolution is required and sufficient computer resources are 
available, a better solution of the master problem is 
necessary, yielding better intermediate lower bounds. 

Besides such a change of input parameters, a structural 
modification was also investigated. Since each capacity 
reallocation routine creates a subproblem solution that is 
feasible and at least as good as the current upper bound, 
this information can be passed to the master problem by using 
this solution as an additional "extreme" point. Experimenta- 
tion on the test problem showed however that no significant 
improvements were obtained. This may be because the dif- 
ference between the interior point solution of the master 
problem and the solution of the capacity reallocation 
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procedure vanishes while the lower bounds are still trailing 
behind. Experimentation on other problems may yield dif- 
ferent results . 

The capacity allocation mechanism can be modified 
further. For RSD (B) , capacity allocation for an interior 
point amounts to a redistribution of the available slack. As 
stated in equation (24) , each commodity receives the same 
additional amount. A proportional allocation would be given 

by ?pj = Xp. + Xpj (b^ - A^^)j / (A^)j. (71) 

The disadvantage of this procedure is that a commodity with 
zero flow on that arc gets zero capacity allocated. To 
overcome this drawback, a convex combination of both methods 
is possible: 

^Pj " ^Pj ^Pj - Ax)j / (A^)j 

+ S>2 (bi - Ax)j / |P| (72) 
where JJ^^, ^2 > 0 and 152+^2”^* 

Experimentation with the test problem 4H showed improve- 
ments in the upper bounds in both cases, but since the lower 
bounds are still weak, the overall gain is not significant 
for this problem. 

RSD(B) would be superior to RSD(P) if we could establish 

A 

better lower bounds. The dual estimates (i are a function of 
the barrier parameter h and the slack on the joint capacity 
arcs. Different dual estimates are obtained if we use a 
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different value of h for each individual joint capacity arc. 
This approach leads to the weighted barrier function like the 
one proposed by Megiddo [Ref. 14] in general linear programm- 
ing : 

B(h,x,w) = cx - h Xj Wj In (b^ - Ax) j , (73) 

where w is any real vector with positive components. Some 
experiments were done with different weights. All arcs that 
are violated in the initial solution with completely relaxed 
joint capacity constraints are more likely to be tight in the 
optimal solution. Therefore they are assigned a reduced 
weight to enable a faster approach to the boundary. Another 
weight factor used was proportional to the remaining slack on 
the corresponding arc. Unfortunately, neither attempt 
provided better solutions . 

The final modification is a hybrid barrier-penalty 
algorithm. Starting with the barrier function decomposition 
in problem lOH, we shift to the penalty algorithm after 16 
iterations using the extreme points generated so far. The 
main problem is the adjustment of the parameter h. Since we 
obtain good lower bounds for relatively small values of h, 
we reset h to the same initial value that we used in the 
independent approach. The results displayed in Figure 13 are 
unimpressive but further experimentation may improve results . 
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Figure 13 : Trajectory of the Hybrid Method 

Finally, we obtain a solution to the 100 commodity 
network flow problem lOOH which is about 3.5% optimal within 
2650 seconds after 25 iterations and about 2.5% optimal 
within 3800 seconds after 38 iterations. The initial value 
h® has been increased to h® = ^^*^max* Staniec [Ref. 4] 
reports a solution to this problem obtained by a penalty 
method that achieved 4 % optimality after 1000 seconds and 
finished with 1.5% optimality after 3000 seconds. Our method 
did not show the same performance. The solution trajectory 
as shown in Figure 14 seems to indicate that the upper bounds 
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are already very tight whereas the lower bounds are again 
poor. The objective function value (not shown) is still less 
than the lower bound, a further decrease of h is necesarry 
for convergence. 




Figure 14 : Solution Trajectory, Problem lOOH 

The distribution of the CPU time between the master 
problem and the subproblems changes as we increase the number 
of commodities. Although the solution of each pure network 
flow problem requires less than one second CPU time per 
comir.odity in the average for all test problems, we find that 
the subproblems are expensive to solve and consume the 
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largest portion of the total CPU time. This distribution is 
shown in Table 2. The amount of time not used by the master 
problem or the subproblem is mostly consumed in the 
input/output operations. 





Problem 4H 


Problem lOH 


Problem lOOH 


Master Problem 


25.5% 


11.2% 


7.1% 


Subproblem 


49.9% 


66.9% 


86.3% 


Other 


24.6% 


21 . 9% 


6.6% 



Table 2; Distribution of CPU Time for Test 
Problems of Increasing Size 

However, the subproblem solves for each commodity indepen- 
dently and more than one pure network flow problem could be 
solved at a time. This feature makes our method highly 
parallel with a potential reduction of nearly 1/|P| in 
elapsed computation time for the subproblem. 
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V. CONCLUSIONS 



The method of logarithmic barrier function decomposition 
is a useful tool for the solution of large-scale multicom- 
modity flow problems. The algorithm RSD (B) is a variant of 
the price-directive decomposition method. It generates a 
sequence of interior points which provide intermediate 
feasible solutions while converging towards the optimum. The 
technique of restricted simplicial decomposition (RSD) proves 
to be useful in the solution of the nonlinear master problem. 
However, RSD (B) seems to be robust and does not require an 
optimal solution to the master problem. It is competitive 
with penalty decomposition methods and relatively easy to 
implement. Since it achieves good feasible solutions at an 



early stage, it 


may 


be used 


as a starting technique in 


other 


algorithms like 


the 


hybrid 


barrier-penalty technique. 


The 


use of RSD(B) 


in 


other 


large-scale multicommodity 


flow 



problems is recommended. 
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